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Abstract 

Quantum Monte Carlo simulations of interacting electrons in solids often 
use Slater-Jastrow trial wave functions with Jastrow factors containing one- 
and two-body terms. In uniform systems the long-range behavior of the two- 
body term may be deduced from the random-phase approximation (RPA) of 
Bohm and Pines. Here we generalize the RPA to nonuniform systems. This 
gives the long-range behavior of the inhomogeneous two-body correlation term 
and provides an accurate analytic expression for the one-body term. It also 
explains why Slater-Jastrow trial wave functions incorporating determinants 
of Hartree-Fock or density-functional orbitals are close to optimal even in 
the presence of an RPA Jastrow factor. After adjusting the inhomogeneous 
RPA Jastrow factor to incorporate the known short-range behavior, we test it 
using variational Monte Carlo calculations. We find that the most important 
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aspect of the two-body term is the short-range behavior due to electron- 
electron scattering, although the long-range behavior described by the RPA 

should become more important at high densities. 
PACS: 71.10.Ca, 71.45.Gm, 02.70.Lq 
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I. INTRODUCTION 



This paper discusses approximate ground-state wave functions for inhomogeneous inter- 
acting many-electron systems such as solids. In particular, we consider wave functions of 
the Slater- Jastrow type, \I/ = e J D, where D is a Slater determinant and J, the Jastrow 
factor, takes account of the electronic correlations. We have two main aims: we wish to de- 
vise a method for generating inhomogeneous Jastrow factors appropriate for use in strongly 
inhomogeneous solids; and we wish to understand the surprisingly accurate results obtained 
when Slater- Jastrow trial functions are used in variational (V) and diffusion (D) quantum 
Monte Carlo (QMC) simulationsSi of weakly-correlated solids such as silicon. Despite the 
apparent simplicity of the Slater- Jastrow form, cohesive energies calculated using VQMC 
are typically an order of magnitude more accuratei~i than cohesive energies obtained us- 
ing Hartree-Fock (HF) or density-functional theory within the local density approximation 
(LDA).i Cohesive energies calculated using DQMC are even more accurate. 

This is not the place for a general introduction to QMC (see Hammond et a/.i for a 
review) , but a brief sketch of the VQMC method may be helpful. The idea is to obtain an 
approximate many-electron ground state by numerically optimizing an explicit parametrized 
trial wave function. Once this has been done, the calculation of expectation values reduces 
to the evaluation of multi-dimensional integrals. Ordinary integration methods using a grid 
become very inefficient when the dimension of the integral is greater than 5 or 10 (equivalent 
to a 2 or 3 electron system in three dimensions), and so the integrals are evaluated using 
Monte Carlo integration.! This approach scales much more favorably with the dimension 
than grid-based integration methods. 

The trial wave functions used in most QMC simulations contain a Slater determinant 
of LDA or HF orbitals and a Jastrow factor that includes pairwise correlation terms, 
Ua utJ j (fi, fj), and one-electron terms, Xo-^fi), where and <ii are the position and spin 
component of electron i. The u and x functions are usually obtained by optimizing spe- 
cific parametrized functional forms according to the variational principle. In simulations of 
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solids, it is common to simplify the optimization by insisting that u be both homogeneous 
and isotropic (i.e., u au(Tj {r^ Vj) is assumed to depend only on the interelectronic distance 

r ij — \ T i ~ r j\)- 

By contrast, our principal aim is to derive a physically- motivated inhomogeneous and 
anisotropic Jastrow factor for nonuniform systems, based on a generalization of the random- 
phase approximation (RPA) of Bohm and Pines. This will enable us to reduce or even 
dispense with the time-consuming optimization procedure. The RPA is known to give 
unphysical results in some cases and is not generally regarded as an accurate quantitative 
method. Here, however, we use the RPA only to guide the construction of an approximate 
trial wave function. Once this wave function has been chosen, energies are calculated using 
the true interacting many-electron Hamiltonian instead of the approximate RPA form. The 
results are therefore variational and much more accurate than standard RPA energies. 

The RPA theory of the homogeneous electron gas is already used in the construction of 
Jastrow factors for QMC simulations.0 This approach predicts that the u function should 
decay like for large r^, but says nothing about the x function (which is zero in a 

homogeneous system) or the short-range behavior of u. Although it is easy to modify the 
u function to make it have the correct cusp-like behavior at short range, the absence of x 
terms implies that the homogeneous RPA Jastrow factor produces inaccurate densities, and 
hence poor energies, when used in strongly inhomogeneous systems. This problem is usually 
fixed by adding a parametrized x function and optimizing the parameters numerically.il 

Here, we generalize the RPA theory to inhomogeneous systems (see Malatesta0 and 
Fahy§ for an alternative approach). We obtain a correlation term that is truly inhomo- 
geneous and anisotropic. Furthermore, as originally noted by Malatesta et ai, 111 we find 
that the inhomogeneous RPA theory automatically produces a Jastrow factor containing x 
terms. The short-range behavior of the inhomogeneous RPA u function is no better than 
that of the homogeneous version but is more difficult to correct. We impose the required 
short-range cusps using a /c-space method which, although approximate, works well. 

An interesting aspect of the inhomogeneous RPA theory is that it provides a justification 
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for using Slater determinants consisting of LDA or HF orbitals. This is common practice 
but is not obviously correct: one might guess that the HF orbitals that are optimal in the 
absence of a Jastrow factor would no longer be accurate in its presence. In fact, the RPA 
theory shows that HF or LDA orbitals are close to optimal whether or not an RPA Jastrow 
factor is present. 

We test our inhomogeneous RPA Jastrow factor by carrying out VQMC calculations for 
a strongly inhomogeneous electron gas. We find that the inhomogeneous RPA \ functions 
are of such high quality that there is no need to resort to the standard but costly numerical 
optimization methods. Surprisingly, however, we gain almost no advantage by using inho- 
mogeneous u functions: a Jastrow factor consisting of a homogeneous u function plus any 
good x function gives almost the same results. Inhomogeneous u functions are often used 
in full core atomic and molecular calculationsBJll but do not seem necessary in the strongly 
inhomogeneous electron gases studied here. Some recent VQMC simulations by HoodEllll 
and Nekovee^l suggest a possible explanation for this observation. Although our system has 
a strongly inhomogeneous exchange-correlation hole n xc (r,r') = [g(r, r') — l]n(r'), the LDA 
seems to provide a better description of the pair-correlation function g(r, r'). The use of a 
strongly inhomogeneous u function may therefore be unnecessary. 

In agreement with previous work@ we conclude that the most important features of 
the Jastrow factor are (i) that the corresponding Slater- Jastrow wave function produces an 
accurate electron density, and (ii) that the u function has the correct cusp-like behavior 
whenever two electrons approach each other. As long as both these conditions are satisfied, 
the calculated energies are normally quite accurate. 

In summary, the five main results of our paper are: 

1. We have extended the RPA to inhomogeneous systems. 

2. We have used the inhomogeneous RPA theory to derive expressions for the x function 
and for the long-range behaviour of the fully inhomogeneous u function. 
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3. We have developed a general method for imposing a cusp on any inhomogeneous 
Jastrow factor expressed in fc-space. 

4. The inhomogeneous RPA analysis has enabled us to explain why Slater- Jastrow-type 
wave functions containing LDA or HF orbitals work so well. 

5. We have implemented and tested the inhomogeneous RPA Jastrow factor and found 
that the calculated one-body term works perfectly, but that the inhomogeneity of the 
two-body term contributes little to the energy at typical valence electron densities. 

The rest of this paper is organized as follows. In Sec. || we describe the Slater- Jastrow 
trial wave functions used in most QMC simulations of atoms, molecules, and solids. Section 



n| presents the RPA theory of the inhomogeneous electron gas and explains how it leads 
to Slater- Jastrow trial wave functions containing both \ terms and inhomogeneous u terms. 
To supplement the somewhat mathematical presentation in Sec. fTJ, Appendix |A] describes 
the RPA in more physical terms. Section [IV] discusses the results of the VQMC simulations 
we have done to test the inhomogeneous RPA Jastrow factor, and Sec. |VJ concludes. 



II. TRIAL WAVE FUNCTIONS FOR QMC SIMULATIONS 

The aim of this paper is to provide a better physical understanding of the success of the 
Slater- Jastrow trial wave functions used in many QMC calculations of atoms, molecules, 
and weakly correlated solids. A Slater- Jastrow trial function is the product of a totally 
antisymmetric Slater determinant D and a totally symmetric Jastrow factor e J . The Slater 
determinant is often split into two smaller determinants, one for each spin value: 

* = e J £) T D i . (2.1) 

This spoils the antisymmetry of the trial wave function on interchange of electrons of oppo- 
site spin, but does not affect expectation values of spin- independent operators.0 Since two 
smaller determinants are easier to deal with than one big one, it also reduces the numerical 
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complexity of the problem. The orbitals used in and are normally obtained from 
LDA or HF calculations. 

The Slater determinants build in exchange effects but neglect the electronic correlations 
caused by the Coulomb interactions. The most important correlation effects occur when 
pairs of electrons approach each other, and these may be included by choosing a pairwise 
Jastrow factor of the form, 



where u aa '(r, r') = u a > a (r f , r). Because the LDA or HF orbitals in and Di already give 
a reasonably good approximation to the density, the introduction of a two-body u function 
usually causes the density of the many-electron wave function to deteriorate. As a result 
the trial energy deteriorates too. It is therefore necessary to introduce one-body x terms to 
adjust the Jastrow factor: 



Note the sign conventions here: our definition of u is the negative of that used by many 
other authors. Most authors also omit the diagonal i = j terms in the sum over % and j, 
but we find it mathematically convenient to include them in our analysis. In homogeneous 
systems the diagonal terms only affect the normalization of the trial function, while in 
inhomogeneous systems they add one-body contributions that may be accounted for by a 
simple redefinition of Xo-( r )- 

Since we are using periodic boundary conditions, it will often prove convenient to express 
the electron density and Jastrow factor in reciprocal space: 




(2.2) 




(2.3) 




(2.4) 



(2.5) 



(2.6) 
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Note that we are using a symmetric definition of the Fourier transformation; the correspond- 
ing back transformation is /(r) = ^ _1/ ' 2 Ek/(k)e" tr k - 

Several different types of Jastrow factor are in common use: 

1. Jastrow factors based on the RPA theory of the uniform electron gas in the form due 
to Bohm and Pines.0 (In this paper we show how to generalization this approach to 
obtain Jastrow factors for strongly inhomogeneous systems.) 

2. Jastrow factors based on the RPA theory of the uniform electron gas in the form 
due to Gaskell.0 This type of Jastrow factor has been widely used by Ceperley and 
coworkers.0 

3. Jastrow factors obtained by the numerical optimization of a parametrized functional 
form. In most cases the optimization is carried out using the variance minimization 
technique in the form developed by Umrigar.0'0 

4. Jastrow factors derived using the Fermi hypernetted chain approximation as given by 
Krotscheck et al@ 



A. The cusp conditions 

In this section we deduce the short-range behavior of u. When the distance = |rj — rj | 
between two electrons approaches zero, the potential energy in the Hamiltonian operator 
diverges like (Except where otherwise stated, we use Hartree atomic units, % = e = 

47T£o — Tn e — 1, throughout this paper.) Since, for any exact eigenstate, Hty is proportional 
to \&, this Coulomb divergence must be cancelled by an equal and opposite divergence in 
the kinetic energy terms. When choosing a trial wave function for use in a QMC simulation, 
it is important to ensure that this exact cancellation still occurs. If it does not, the local 
energy El = H^f/^f, which is the quantity actually sampled in the simulation, diverges as 
r« approaches zero, causing numerical difficulties. 
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Since, away from the nuclei, the single-particle orbitals are smooth functions of the 
electronic coordinates, the determinantal part of the trial wave function depends smoothly 
on the positions of the electrons. For small r^, it may therefore be expanded in the form: 

(D^D l )\ ri ^ Q ^c + a -v ij + ... . (2.7) 

The application of the kinetic energy operator to this series produces a smooth function of 
Tij, whereas the application of the potential energy operator places a divergent factor 
in front of every term. 

It is clear that a trial function containing only the Slater determinants does not show 
the required cancellation of potential and kinetic energy divergences. The cancellation may 
however be imposed by modifying the terms of the expansion as follows: 

c — > (1 + \rij) c , 

(2.8) 



a • r 



-> (1 + ^Tij) a • 



Each term in the series is multiplied by a factor, 1 + ar^, which includes a cusp at = 0. 
The numerical constant a in front of the cusp is 1/(2 + 2n), where n is the order of in 
the original expansion. 

A convenient way to introduce a cusp-like term is to include it in the Jastrow factor e J 
in such a way that the cusp terms (1 + \rij) or (1 + \rij) appear in the expansion of e J 
for small r^. As we only have one Jastrow factor we can only hope to deal with one of the 
terms in the expansion of D^D^ correctly; we choose the lowest order term that is nonzero, 
as this is the one that causes the divergence in the local energy E L . If electrons i and j 
have antiparallel spins, the lowest order term is the constant c; for parallel spins, D^D^ is 
an antisymmetric function of and so the lowest order term is a • r^. The required cusp 
may thus be introduced by imposing the following conditions on n ffi(T) : 



du U (Ti,Tj) 

for antiparallel spins, and 



drij 



\ (2-9) 
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du 



dry 



(2.10) 



for parallel spins.@H 



B. The homogeneous RPA correlation term 

As will be explained in Sec. |T1|, the RPA theory of Bohm and Pine suggests that 
the long-range behavior of the correlation term u in a homogeneous electron gas of number 
density n should take the form: 

u aitTj {r h Tj) = u . a . (r -) = — — , (2.11) 

p y 



where u p = y/Ann is the plasma frequency. This spin-independent two-body term has no 
cusp and is only expected to be correct for large r^. Multiplying Eq. (|2.11|) by 1 — e~ ry / i?0 '* j 
yields 



— (l - e-^'^i) , (2.12) 



which approaches Eq. (|2.11|) for large ry and also has the correct cusp behavior if F ai(Jj 
is chosen appropriately. Because of the spin dependence of the cusp conditions, we need 
different values of F for parallel- and antiparallel-spin electrons. By expanding Eq. ( |2.12| ) 
to first order in r^, we see that the cusp conditions Eq. ( p.9[ ) and Eq. ( |2.10| ) become 

1 



1 _ du u 



for antiparallel spins, and 



2 dr 



1 diijj 



r=0 2F^LU p 



(2.13) 



dr 



(2.14) 



for parallel spins. 
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C. The x function 



The x function needed to counteract the density-adjusting effects of the u function is 
normally obtained by numerical optimization. If we wish to avoid this costly procedure we 
have several options: 

1. We could use no one-body term and hope for the best. This is the correct thing to 
do in a homogeneous system and also proves satisfactory in cases when the spatial 
variation of the charge density is much more rapid than that of the u function. 



2. We could use the result of Sec. |TJ and include a one-body term of the form: 



XM = - E <w(k,k')rv(k') . (2.15) 

k',cr' 

3. We could use a one-body term as given by Malatesta et a/.@ 

Xa(k) = -V^$>-'( k K'( k ) • (2-16) 

cr' 

In cases when the u function is homogeneous, 

u<j<j< (r, r') = *w (r - r') , (2.17) 

we get 

W(k,k') = VF?w(k)5 k , k , (2.18) 



and so Eq. ( |2.15| ) reduces to Eq. ( |2.16| ). Option (2) therefore reduces to option (3). Strictly 



speaking, however, u is never exactly homogeneous unless the electron density is exactly 
uniform, in which case both options (2) and (3) merely state that \ = 0. 
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III. THE RANDOM-PHASE APPROXIMATION FOR INHOMOGENEOUS 

SYSTEMS 



This section starts with a brief introduction to the ideas behind the standard RPA treat- 
ment of homogeneous systems.0110 We then generalize the RPA theory to inhomogeneous 
systems. This generalization shows us how to construct approximate ground-state wave 
functions of the Slater- Jastrow type. The RPA forms of the u and x functions are closely 
related and both depend on the electron density. 

A. Review of the RPA for homogeneous systems 

Let us first give a brief overview of the standard RPA as formulated for homogeneous 
systems. The starting point is the usual Hamiltonian, 

H = - 2 T,V l +2^— ^E^. (3-1) 

of a uniform electron gas in a volume V. Note that we have adopted a slightly unusual 
definition of the number density operator, 

^k = J=£^, (3.2) 
Vv i 

which includes a l/y/V factor in order to be consistent with the symmetric definition of 
the Fourier transform used throughout this paper. The /c-points are chosen in accordance 
with the geometry of the system, which is taken to obey periodic boundary conditions. The 
system is assumed to be charge neutral and so the k = terms are omitted from the fc-space 
summations. 

As we are interested in describing the long-range correlations due to the electron-electron 
interaction, we wish to make a connection between the Hamiltonian of Eq. (|3.1|) and the 
long wavelength density fluctuations known as plasmons. Our eventual aim is to split the 
Hamiltonian into an electronic term with short-range interactions and a plasmon term that 
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is only weakly coupled to the electrons. Once the electron and plasmon parts of the Hamilto- 
nian have been (almost) decoupled in this way, we shall see that the long-range correlations 
are described by the ground-state wave function of the plasmon part, which is simply a 
set of quantum mechanical harmonic oscillators. At shorter wavelengths the plasmons are 
not well defined and the collective plasmon description of electronic correlation ceases to be 
valid. Instead, the electrons feel a short-range screened interaction that produces additional 
electron-electron scattering-like correlations. 

The assumption of weak plasmon-electron interaction is reasonable at small k since long 
wavelength plasmons are long lived. For larger k values, however, the almost flat plasmon 
dispersion curve runs into the continuum of electron-hole pair excitations and the plasmons 
are no longer well defined. In a uniform electron gas this happens at a wave vector k c given 
by§ 

k c « X -k F r x J 2 oc n 1 ' 6 , (3.3) 

where kp is the Fermi wave vector, r s ao is the radius of a sphere containing one electron on 
average, and ao, the Bohr radius, is the atomic unit of length. This estimate of the cutoff 
should also be applicable to inhomogeneous systems as long as the density does not vary 
too much. We see that for typical metals with r s values of 2 or 3 the cutoff is of the order 
of the Fermi wave vector. We also see that the higher the density the better the plasmon 
description should be. The inverse of k c is a measure of the minimum length scale on which 
the electronic correlations may be described in terms of plasmons. 

The first step in the derivation of the homogeneous RPA0 is to introduce a new pair 
of conjugate operators, 7r k and g k , for every wave vector k with modulus k < k c . These 
operators act in a new Hilbert space that we call the oscillator space, and transform like 
under Hermitian conjugation: 7T k = 7rl k and g k = gl k . The physical and oscillator Hilbert 
spaces are quite distinct, and so 7r k and g k commute with fj and p>j. For the time being 7r k 
and g k have little physical meaning, but later in the derivation they will become associated 
with the plasmon coordinates. 
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The oscillator-space operators are now used to define a new Hamiltonian, 
^ BP = 2? P ^ +27r ?^ \T^¥ 

+2 E ^4 - E (^) ^4 , (3-4) 

which acts in an enlarged Hilbert space that is the product of the physical space and the 
oscillator space. The Bohm-Pines Hamiltonian H^p may be written in the form 




and so its eigenvalue spectrum is bounded below. 

We now restrict our attention to those states |$) from the enlarged Hilbert space that 
satisfy the subsidiary condition: 

7r k |$) = for all k < k c . (3.6) 

Any state obeying Eq. ( [3.6D may be written in the form |<&) = |?/>)|7r=0), where is a 
state in the physical Hilbert space and |tt=0) is the oscillator-space state satisfying 

7r k |7r=0) = for all k < k c . (3.7) 

(The bold symbol w is shorthand for the vector of all the tt^ with k < k c .) The set of 
product-space states satisfying the subsidiary condition may therefore be put into one-to- 
one correspondence with the set of physical states. Equally, given any eigenstate with 
eigenvalue E of the original Hamiltonian, the product state |$) = \ip)\n=Q) satisfies the 
subsidiary condition, Eq. (|3.6| ), and is an eigenstate of Eq. ( |3.4j ) with the same energy. This 
is because all the new terms in involve 7r k and hence give zero when operating on |7r=0). 
As long as the subsidiary condition is obeyed, the additional degrees of freedom may simply 
be regarded as dummy variables, both in the wave function and the Hamiltonian. 

It will be proved later on in this paper that the overall ground state |$o) of the ex- 
tended Hamiltonian, Eq. ( |3.4|) , automatically satisfies the subsidiary condition, Eq. ( p.6|) . 
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The physical-space part of the ground state of Hbp is therefore the ground state of the orig- 
inal Hamiltonian, and we might as well study the ground-state properties of the extended 
Hamiltonian as those of the original. 

The next stage in the derivation is to make a unitary transformation, which will be de- 
scribed in more detail in Sec. [Ill C| . Once the extended Bohm-Pines Hamiltonian has been 
transformed and several supposedly small terms have been dropped, the transformed Hamil- 
tonian may be written as the sum of a spatial Hamiltonian with a short-range interaction, 

fj 1 , o„ \- Kfj 2vriV_ 1 

H sr = 7>2^Pi + 2?r —2 ^^-p' (3 ' 8) 

and an oscillator-like plasmon term, 

Hp = 2 ^ (^k + ^QkOk) > (3-9) 

where lu p = (AtcN/V) 1 ^ 2 is the plasma frequency. One of the approximations made during 
the derivation of these results involves the replacement of an exponential factor of the form 
exp[i(k — k') ■ by its average value. Since the electronic positions r$ are effectively random 
in a homogeneous system, the phase of the exponential is also random unless k = k', and 
so the average is 1 if k = k' or zero otherwise. This is the eponymous random-phase 
approximation. 

The approximate Hamiltonian, 

Hrpa = H sr + Hp , (3.10) 

still includes short-range electron-electron interactions and so cannot be solved exactly. We 
may, however, use the HF method or the LDA to obtain approximate eigenstates of H sr . 
These approximate eigenstates, which are of course Slater determinants, may then be mul- 
tiplied by eigenstates of the oscillator Hamiltonian, Eq. Q3.9Q , to obtain approximate eigen- 
states of the full Hamiltonian -^rpa- The unitary transformation may then be inverted to 
obtain approximate eigenstates of i^BP- Unlike the exact ground state, the resulting approx- 
imate ground state of i^BP does not satisfy the subsidiary condition, Eq. (|3.6|) , exactly. It 
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may, however, be projected onto the subspace of states that do satisfy the subsidiary condi- 
tion to obtain an approximate ground state in the physical Hilbert space. This approximate 
physical ground state has the homogeneous Slater- Jastrow form, 



V>o({xi}) = exp 



2tt v n k nj 

U P k<k c K 



D{{^}) , (3.11) 



where = V~ 1 l 2 Yl,i^' Vi and Xj = (r^, cr,) describes the position and spin component of 
electron i. The RPA u function is therefore independent of spin and proportional to l/r^ 
at large r^. As noted earlier, the cutoff wave vector k c is usually chosen to be of the order 
of the Fermi wave vector kp. 

B. The RPA Hamiltonian in inhomogeneous systems 

The aim is to find an approximate ground state |?/>o) with energy E of the following 
Hamiltonian: 

# = ^Ep 2 + E^) 

, 9 ^y-%4 2ttA_ 1 

where V(r) is an applied potential. The "physical" Hilbert space in which this Hamiltonian 
acts will be denoted Hr. For reasons that will become clear later on, it is sometimes 
convenient to rewrite H in the form, 

# = ^Ep 2 + E^) 



, o K™1 2ttA _ 1 



(3.13) 

where V(r) is defined via: 
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V(r) = V(r)-^--AV(r), (3.14) 

ab^E^A- ( 3 - 15 ) 

AV,(r)= -£(S) (3 ' 16) 

and the 7r k are arbitrary numbers satisfying 7r k * = 7r° k . 

Let us now introduce conjugate pairs of operators, and q^, acting in a different Hilbert 
space Ho, which we call the oscillator space. The "momentum" operator 7r k and the "posi- 
tion" operator g k are taken to be the Fourier transforms of field operators 7r(r) = 7r^(r) and 
q(r) = gt( r ) : 

& = -j=J v e*~q1(r)<Pr, (3.17) 
n k = ^= I e ik - r 7r(r)d 3 r , (3.18) 



which satisfy the commutation relation 

[7r(r),g(r / )] = -^(r-r'). (3.19) 

It therefore follows that 7r k = 7r_ k and g k = g_ k obey 

[TTk, <7k'] = -^k,k' • (3.20) 

We can form an extended Hilbert space by taking the product space 7i cx t = Hr®Ho- 
If we denote the identity operators in the real and oscillator spaces by 1r and lo respec- 
tively, we can define extended-space operators such as f cxt = tr <g> 1q and g ex t = Ir <8> <?o 
corresponding to any operator belonging to one or other of the constituent Hilbert spaces. 
It is obvious that all the extended "i?" operators commute with all the extended "O" oper- 
ators. From now on we shall omit the identity operators and drop all "O" , "i?" , and "ext" 
subscripts. Operators will always be denoted using hats, and so anything without a hat may 
be assumed to be a (possibly complex) number. 
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Now consider the Hermitian extended Hamiltonian, 



i i 

k K k K 



+ 5E^k- E fS) 1/2 ^4, (3-21) 



2 ^ ^ y^2 



which is an operator in the extended Hilbert space 7Y ex t- Note that the potential is still 
the same as given by Eq. (|3.14|) ; the variables 7r£ that were used to construct V have not 



been replaced by operators in Tio-, but remain ordinary complex numbers. Because Hbp 
and 7r commute they may be diagonalized simultaneously, and hence all eigenstates of the 
extended Hamiltonian H^p may be written in the form \i/j^)\tt), where 7Tk|7r) = 7rk|7r) for 
all k < k c . 

Let j^o) and E be the ground-state wave function and ground-state eigenvalue of the 
physical Hamiltonian H. If we define an eigenstate |7r°) of the fc^ operators such that 
^kl 71 " ) = ^kl 71-0 ) f° r a ^ & < ^ci it follows that \ipo)\7r°) is an eigenstate of H^p with the same 
eigenvalue E . It need not, however, be the ground state of Hpp, which might correspond 
to a different wave function I^Vm) 1 7r mm ). The question that now arises is how to choose the 
constants 7r£ (and hence the modified potential V) such that \%1)o)\it ) is in fact the ground 
state of f^BP and not just some other eigenstate. In other words, we have to find a link 
between the 7r£ and |^o)- This can be achieved using the Hellmann-Feynman theorem. 

Consider the lowest energy state !$„-) = \ip n )\ir) corresponding to some fixed oscillator- 
space eigenstate \ir). The energy eigenvalue of will be denoted E*. When H BP acts 
on |$7r), the operators 7Tk m ay be replaced by their eigenvalues 7r k , and hence j^) is in fact 
the ground state of 

^ = ^Ep? + E^) 

i i 

\~ 2 ^\k~ 2 
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1 /4tt\ 1/2 
+ g E - £ ( p ) Tk4 • (3-22) 

Note that operates not in the extended Hilbert space but in T~Lr. Furthermore, it is 
important to keep in mind that V is formed using the as yet unknown numbers 7r£. The 
overall ground state of H-gp corresponds to the minimum of with respect to 7r. As we are 
now looking at a standard quantum mechanical problem formulated in the physical Hilbert 
space Hr, the Hellmann-Feynman theorem gives: 

dH 7 



7T-k 



^ (3.23) 



J <tM4ltM , (3-24) 



mm 

> 



where is assumed normalized. If the overall ground state of if bp occurs when w = 7r 
it follows that: 

= *T ~ (p ) (Vv^i4l^> • (3-25) 

We wish to choose the constants 7r£ such that 7r£ = vr^ 1111 , since in this case we have 
already argued that the physical-space part l^o) of the Bohm-Pines ground state 

|$ 7r o) = |^ 7r o)|7r ) (3.26) 

is equal to the exact physical ground state |^o)- The sought after link between |^ ) an d the 



7r£ is therefore: 



/4vr\ 1/2 
k 2 



n- k )o , (3-27) 



where (n k )o = (V'ol^klV'o) is a Fourier component of the ground-state electron density n(r). 
Note that Eq. (|3.27 ) is not an operator equation; it simply links one number, 7r°_ k , to another, 



(^-k)o- Since all the #k operators have k < k c , Eq. ( |3.27| ) only applies when k < k c . 

We still have to verify that the stationary point 7r mm is in fact the minimum of E 7 *. 
Eq. ( p.24| ) determines the slope at any given 7r. The matrix of second derivatives of E 77 , 



which determines the curvature, may therefore be obtained by differentiating Eq. ( ft.24 ): 
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= far 1 (-k,-k')/A;', (3.29) 



where 



x(k,kO = -4=l^ (3-30) 

V 47T C7T k / 



is the static susceptibility matrix and e -1 (k, k') is the inverse dielectric matrix. If we make 
the reasonable assumption that all the eigenvalues of the dielectric matrix are positive,@~i2l 
it follows that all the eigenvalues of the matrix of second derivatives of E n are also positive. 
The stationary point given by Eq. fl3.25| ) is then a minimum. 

Eq. ( 3.25|) is the generalization to inhomogeneous systems of the "subsidiary condition" 



of Bohm and Pines.0 It can be rewritten as 

fi k |$) = (k < k c ), (3.31) 

where 

= vr k - [-) (n k )ol , (3.32) 

and has to be obeyed by |<&o)> the exact ground state of Hbp- In the case of a homogeneous 
system, Eq. ( |3.31| ) reduces to 7r k |^) = as derived by Bohm and Pines.0 

The effect of the subsidiary condition is to reduce the number of degrees of freedom of the 
extended Hilbert space TL ex t- The reduced Hilbert space satisfying the subsidiary condition 
is spanned by the set of eigenstates of Hbp of the form \if))\7r), where 

TTk = ( tj J (n k )o for all k < k c . (3.33) 



If we have set the parameters 7r k using Eq. (|Q7D , the states are then eigenfunctions of 



the original Hamiltonian, Eq. ( p. 12] ). This follows from the definition of the potential V"(r), 
which is such that the Hamiltonian in Eq. (|3.22|) is the same as the original Hamiltonian of 



Eq. ( |3.12|) when the value of it is consistent with the subsidiary condition: H n | o — H. 
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The subspace singled out by the subsidiary condition Eq. (|3.31| ) is thus equivalent to the 
original Hilbert space 7Yr. 

Let us now use the subsidiary condition to evaluate AE and AV(r). For AE we get 

= -J n l {r)n\r')W{r - r')d 3 rd 3 r' , (3.34) 

where n l (r) is the long wavelength (k < k c ) part of the ground-state electron density and 
W is the periodic (Ewald-summed) Coulomb interaction. The constant AE is therefore the 
long wavelength contribution to the Hartree energy. For AV we get 

A^r) = - J= E ^Oo^ , (3-35) 

V * k<k c 

and hence 

V 2 [AV(r)] = 4™'(r) . (3.36) 

AV(r) is therefore the Hartree potential corresponding to the long wavelength Fourier com- 
ponents of the electronic charge density. Because AV(r) is subtracted from V(r) to give 
V(r), the extended Hamiltonian Hbp contains a reduced external potential. The long-range 
part of the mutual repulsion of the electrons has been absorbed into the plasmon degrees of 
freedom via the subsidiary condition. 

C. The unitary transformation 

We have now concluded that we can concentrate on the ground state of the Bohm-Pines 
Hamiltonian, Hbp, from Eq. ( |3.21|) instead of the ground state of the original Hamiltonian 



from Eq. (|3.12|) , provided we choose the constants 7r£ and hence the effective potential V 
in accordance with Eq. ( |3.27| ). The oscillator-space part of the ground state |$ ) of H B p is 
then equal to \tt°), and the real-space part is the physical ground state \ipo)- 

Eq. ( |3.27| ) specifies 7r£ in terms of the ground-state density, which is obtained by solv- 
ing if bp- Unfortunately, this Hamiltonian depends on the parameters 7r£ calculated from its 
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ground state, and so we are faced with a self-consistency problem analogous to those encoun- 
tered in bandstructure calculations. We could in principle devise an iterative algorithm to 
home in on a self-consistent solution, but this is unnecessary. We only need the ground-state 
density, not the wave function itself, and it is known that the LDA gives reasonably good 
ground-state densities in most solids. In practice, therefore, we can use the LDA density 
(^k)o DA to obtain a good approximation for 7r£: 

-47ry/ 2 



If. \ LDA 
k 2 ) ™° 



(3.37) 



Unsurprisingly, the Bohm-Pines extended Hamiltonian cannot be solved exactly. It may, 
however, be solved approximately by means of a unitary transformation.0 We use the unitary 
operator 



S = exp 



- % 2. ( T2 ) ^ llk 

k<k c 



(3.38) 



to transform an eigenstate |$) of H^p into an eigenstate 

|$ ncw ) =S\$) 



(3.39) 



of the transformed Hamiltonian H^T = SHbpS'. (In general, all operators transform 
according to new = SOS^.) The position operators f j and g k are unchanged by the trans- 
formation because they commute with S. The momentum operators transform as follows: 



Pi 



Pi 



" s "new 



Pi + % 



7T k 



ik-Ti 



k 2 ) 



n k , 



(3.40) 
(3.41) 



where £ k = k/k is a unit vector in the k direction. The transformations of the momentum 
operators may be checked using the general expansion, 



e x Oe 



-x 



+ 



X,0 



X, 



X,0 



+ ... 



(3.42) 



with S = e x . When O is one of the momentum operators, 7r k or pj, the first commutator 
[X, O] only contains position operators, and so higher order commutators such as [X, [X, O}} 
all vanish. 
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ik f,; 



The final result of the unitary transformation defined by Eq. ( |3.38|) is the Hamiltonian: 

H bp = 2^Pi + 27r -p- 
2vriV 1 ~ 1>r , 

+i (^) 1/2 £? £k '(^"9 &e 

+ -/= E ( £k ' £k ') <?k<?-k'™k-k' 

v K k,k'<k c 

+ ^ E ^7T-k , (3.43) 

which is obtained by replacing ]3j and 7i"k in Eq. fl3.21|) by pf ew and #£ ew . The subsidiary 
condition becomes 

^ ew |$ ncw ) =0 {k<k c ) , (3.44) 



ni cw = 7r k + f ^ V2 (n k - (n k )o) • (3.45) 



where 



\k 

We now make the random-phase approximation, which amounts to replacing the nk_k' = 
V^ 1 / 2 J2i e^ k_k ')' ri factor in the fourth line of Eq. ( |3.43| ) by its ground-state expectation value. 
In uniform systems the electronic positions are random and so the phases are also random; 
the expectation value of rik-k' is therefore equal to iV5k,k' / VV. In inhomogeneous systems 
we have to evaluate the expectation value of the (untransformed) operator nk-k' i n the 
transformed ground state |$ new ). Since the density operator n(r) = ~ r ) commutes 

with the unitary transformation, it follows that 

($ ncw |n(r)|$ new ) = ($|£tfc( r )£|$) = ($|n(r)|$) . (3.46) 

The required expectation value of rhc-k' is therefore equal to the Fourier component (nk-k')o 
of the ground-state electron density of the original (untransformed) Bohm-Pines Hamiltonian 
from Eq. ([[21]). 
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The term on the third line of Eq. Q3.43T ) may be rewritten as 



-) 2^^k^k-Jk, (3.47) 



2 \ V J 



where j k = J2i{~Piy e lk ri } is the current density operator. Following Bohm and Pines,H this 
plasmon-electron coupling term will be neglected. To justify this approximation (and indeed 
the RPA) we can appeal to the measured physical properties of interacting electron gases; we 
know that the plasmons are well defined when k < k c , and hence that the plasmon-electron 
coupling terms must indeed be small. 

A more physical discussion of the RPA may be found in Appendix 

D. The RPA ground state 

The two approximations described above decouple the electrons and plasmons and reduce 
the transformed Hamiltonian of Eq. ( p. 43 ) to the RPA Hamiltonian Hrpa = H sr + H p . The 



first two lines of Eq. ( 3.43|) yield the short-range electronic Hamiltonian, 



2?riV E4 + E^), (3-48) 



V ^ k 2 

and the last two lines yield the plasmon Hamiltonian, 

H p = -(TT-Trt + q-M-qt) , (3.49) 

where the matrix M is given by 

M kjk , = (e k • Sk ,)i / e l ^'> r ul(r)d 3 r , (3.50) 

and we have introduced a position dependent local plasma frequency defined by 0Jp{v) = 
47m(r). The full ground state of Hrpa is the product of the ground states of H sr and H p . 
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1. The plasmon ground state 



If we choose to work in a representation in which the 7Tk operators are diagonal, the 
plasmon ground state takes the standard simple harmonic oscillator form: 



fy„ oc exp 



I E K(m-^ 



k,k' 



7T k / 



(3.51) 



The matrix M" 1 / 2 is well defined since all the eigenvalues of the Hermitian matrix M are 
greater than zero. The fact that the 7Tk and operators are non-Hermitian may cause some 
confusion here, but one can easily rewrite the /c-space Hamiltonian of Eq. ( |3.49|) in real 
space using Eqs. ( |3.17|) and (|3.18|) . The real-space operators 7r(r) and g(r) are Hermitian, 
and so the plasmon Hamiltonian is then a set of coupled Hermitian harmonic oscillators. 
The ground state of these oscillators is 



oc exp 



/ 7T 

2 



(r)M- 1/2 (r,r , )7r(r , )ci 3 rdV 



(3.52) 



which reduces to Eq. (|3.51|) when re-expressed in fc-space. 



2. The short-range ground state 
If we make use of the expression for V from Eq. ( |3.14j ) and the condition 7r^ 



47r//c 2 (nk)o from Eq. ( |3.27| ), the short-range Hamiltonian of Eq. fl3.48| ) becomes 



^ = ^Ep*+E^) 

z i i 
, o \- %4 27riV _ 1 



+ 47r(n k ) 



-ik fi 



- 2n E ^>°^>° . (3.53) 

The first two lines are identical to the original Hamiltonian, Eq. (|3.12| ), but with the small 
k (long wavelength) contributions to the electron-electron interactions omitted. The third 
line is the Hartree potential corresponding to those long wavelength Coulomb interactions, 
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Vv 



1 



EE 



% A/ fee 



47r(r2 k ) 
fc 2 




(3.54) 



and the fourth line is the Hartree energy, which is subtracted to prevent double counting. 
The short-range Hamiltonian is therefore equivalent to the original Hamiltonian, 



the Hartree approximation. Since the long wavelength parts of the effective potentials used 
in Hartree-Fock and LDA calculations are both dominated by the Hartree contributions, 
we might equally well say that H sr is equivalent to the original Hamiltonian but with the 
small k Coulomb interactions approximated using Hartree-Fock or LDA, provided the HF 
or LDA densities are sufficiently similar to the exact ground state density. The short-range 
Hamiltonian still contains the full Coulomb interaction for k > k c , and so still diverges like 
1/Vjj whenever two electrons approach each other. The electron-electron cusps therefore 
appear in the short-range electronic wave function, not in the Jastrow factor that describes 
the plasmons. 

In practice, of course, we do not attempt to solve H sr exactly, but treat it within an 
independent electron approach such as Hartree-Fock or LDA. This additional approximation 
replaces the short-range part of the electron-electron interaction by a mean field, which 
simply adds to the long wavelength mean field already introduced by the RPA. The overall 
effect is equivalent to starting from the original Hamiltonian and replacing the full interaction 
by a mean field. This implies that one can obtain the short-range "electronic" part of the 
RPA wave function by starting from the original fully interacting Hamiltonian and treating 
it using any sensible mean-field approximation. The best single-particle orbitals to use in 
the Slater determinant are therefore very close to the familiar Hartree-Fock or LDA orbitals; 
they are not significantly altered by the presence of the RPA Jastrow factor from Eq. ( |3.61| ). 

One drawback of treating the short-range Hamiltonian within a mean-field approximation 
is that this neglects the electron-electron cusps that should be present in the short-range 
electronic wave function. The cusps play an important role in reducing the total energy of 
the many-electron system, and so the trial wave function may be significantly improved by 



Eq. ( |3.12j ), but with the long wavelength parts of the Coulomb interaction treated within 
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building them into the Jastrow factor. 



E. Inverting the unitary transformation 

The ground state of Hrpa is the product of the ground states of H sr and H p , neither of 
which commutes with the transformed subsidiary condition, Eq. ( |3.45| ). This implies that, 
unlike the ground state of H&p w , the ground state of H-rpa need not obey the subsidiary 
condition automatically. In consequence, the approximate ground state of Hpp obtained by 
applying the back transformation, S\ to the ground state of Hrpa, need not be an eigenfunc- 
tion of the plasmon momentum operators, and we can no longer extract an approximation to 
the spatial ground state by simply forgetting about the \ir) factor in a product wave function 
of the form \i/j)\tt). Fortunately, however, the subsidiary condition is still exact (no approx- 
imations were made in transforming it), and so still defines the subspace of the extended 
Hilbert space in which the true ground state lies. We can therefore take the ground state 
of the approximate Hamiltonian, Hrpa, and project it onto that subspace. The projection 
operator may be applied before or after the back transformation, but if we choose to make 
the back transformation first it is not difficult to see that the required projection operator 



is 



II h< = 7r°)<7r k = 7r°| . (3.55) 

As discussed in Sec. [Ti l D 2| , we approximate the ground state of H sr as a Slater determi- 
nant D, and so the approximate ground state of the full Hamiltonian Hpj>a is $rpa oc ^ p D. 
We can now obtain an approximate ground state of the original Hamiltonian ifpp by back 
transforming using the inverse of the unitary transformation. The only important effect of 
the back transformation is to shift the numbers 7r k appearing in by — (An /k 2 ) 1 ^ 2 ^: 



7T k -> vr k - [An/k J n k , (3.56) 
where n k = V~ 1 / 2 J2i e ik ' r *- This can be verified by observing that, when evaluating a back- 
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transformed wave function \l/ old ({rj}, {7r k }) = {7Ti c }|S , '''|\l/), we can apply the transfor- 

mation to the bra ({rj}, {7r k }| rather than the ket \^f). But since 

TtkS \{ri}, {7r k }) = SS^kSUri}, 

S(n k -(4n/k 2 ) 1/2 h^\{r t },{n k }) 

: /2 y 

we see that 



7T k - (A7f/k 2 ) n^SUn},^}) 



S\{ri},M) 



{r 4 },{7r k -(4vrA 2 ) 1/2 n k }) . 



(3.58) 



1 /2 

The 7Tk eigenvalues of the transformed bra are therefore shifted by — (Air/k 2 ) n k relative 
to those of the orig inal bra. As a result, ^ old ({ri}, {7r k }) = ^({r,}, {7r k - (An/k 2 ) 1/2 n k }). 

Applying the projection operator given in Eq. (|3.55|) replaces the remaining 7r k by 7r£ = 
(47r/A; 2 ) 1 / 2 (ri k )o. (In the homogeneous case this is zero.) All in all, then, the spatial part of 
the approximation to the ground state is 



* oc $jD = exp 



D 



where 



u(r, r') = —An 







[m~¥ 






k,k' 



Ak'-r' 



V 



N 



kk' 



N 



(3.59) 



(3.60) 



The RPA Jastrow factor includes constant terms, one-electron terms, and two-electron 
terms. The constant terms may be ignored as they only affect the normalization of the wave 
function. The remaining one- and two-electron terms may then be disentangled and the 
Jastrow factor rewritten in the form, 



j oc exp 



i,j * 



(3.61) 



28 



where u(r, r') and x{ r ) are defined via: 



u(r, r 



k,k'<k c 



-ik-r 



47T 



and 



47T 



^ k,k'<k c 



kk' 



M~ l/2 



(3.62) 



k,k' 



fcfc' 



(n k ')o 



- / «(r,rV(r')dV , 



(3.63) 
(3.64) 



where, as in Eq. ( |3.34| ), n (r) is the long wavelength k < k c part of the ground-state electron 
density. The derivation of Eq. ( |3.63| ) made use of the symmetry M^k' = A^_k',-k, which 
follows from Eq. (|3.50|) . In /c-space, the relationship between u and x takes the form, 



x(k) = - E u(k,k')n(k') 



(3.65) 



discussed in Sec. II C 



In a homogeneous system, Eq. (|3.50|) states that M^k' = ^pO~k,k'- The x function therefore 
vanishes and the u function becomes 



u 



horn i 



k,k') 



4tt 1 

Ok k ; 

u n kk' ' 



(3.66) 



Transforming to real space we obtain 



M hom (r,r') = u hom (|r-r'|) 



horn / 



ikr+ik'r' ^ r 

kk' ' 



11 r - 

'vZT E e 

J_J_ c -»k-(r-r')47T 



(3.67) 



If k c is set equal to infinity this gives 

w hom (|r-r'|) 



oo p \r — r 



(3.68) 



For finite k c , the divergence of u(\r — r'\) at small |r — r'| is suppressed, but the l/\r — r'\ 
decay at large |r — r'| remains more or less unaltered. 
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F. The cusp conditions in inhomogeneous systems 



Section 11 B explained how cusps may be built in to a homogeneous RPA Jastrow factor 
by adding an exponential factor to the u function: 

u(r) = — —(l-e- r/F ) . (3.69) 

At large r this u function has the l/(u> p r) behavior implied by the RPA, while at small r it 
tends smoothly towards the required cusp at r = 0. When supplemented by appropriate \ 
functions, such Jastrow factors are remarkably successful. It is therefore worth considering 
how we might add cusps to our inhomogeneous RPA u function. 

This is not easy, since the inhomogeneous u function is given as a complicated truncated 
double Fourier series. The series determines the behavior of u(r, r') when r — r' is large, 
and we have to find a way of splicing this known long-range behavior onto the cusp at 
small r — r'. The cusp fixes the slope of u(r,r') as |r — r'| — > 0, but does not determine its 
position-dependent value at the point r = r'. This makes it difficult to implement simple 
interpolation schemes that use different functions to describe u at small and large r — r'. 
The introduction of a multiplicative factor, as in the homogeneous u function of Eq. ( p.69| ), 
is equally problematic. 

It turns out that this interpolation problem is easiest to handle when expressed in k- 
space. This might seem unlikely at first, since a true cusp can only be generated by a 
computationally intractable infinite Fourier sum. In practice, however, it appears that satis- 
factory approximate cusps can be introduced using Fourier sums with manageable numbers 
of terms. The fc-space cutoff needs to be large enough to ensure that the real space volume 
within which the cusp is not represented correctly is small, and hence that the resulting 
errors contribute little or nothing to expectation values. 

Eq. (|3.69| ) can be Fourier analyzed to give: 

M(k) = ^ P ¥WTUW) ■ < 370 > 
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Eqs. (|2.13|) and ( |2.14|) show that 1/F 2 = Cu p , where C = 1 for antiparallel spins and 
C — 1/2 for parallel spins. Hence 

/i \ -47T C 

n ( k ) = —i= r 0,70 — — r • (3-71) 
W k 2 (k 2 + Cup) v 1 



In section |LV B| we test a homogeneous m function defined using a truncated Fourier series 



of this form, and find that most of the cusp energy can be retrieved using a reasonably low 
cutoff. 

Eq. ( p. 71 ) defines a natural /c-space crossover, k x , given by k 2 = Cu p . The terms with 



k < k x produce the RPA behavior at large r, while for k > k x we have u(k) oc l/k A , which 
generates the cusp. The density dependence of k x differs from that of the plasmon cutoff k c 
from Eq. ( |3.3|) . It turns out, however, that for typical metallic densities k c and k x are both 
of the order oikp- The k 2 + Cuj p factor in the denominator of Eq. (|3.71|) therefore allows us 
to introduce the cusp without significantly affecting the large r (k < k c ) behavior implied by 
the RPA. If the density is extremely small, k x (oc u^ 2 oc n 1 / 4 ) is smaller than k c (oc n 1 / 6 ), 
and so the fc-space method of imposing the cusp is no longer consistent with the RPA limit 
we expect when k < k c . 

Equation (|3.71| ) suggests a simple fc-space prescription for building a cusp into the inho- 
mogeneous Jastrow factor. We write the inhomogeneous Jastrow factor as a double Fourier 



Mir, r 



' -^£e- lk - r M(k,kV k '- r ', (3.72) 



noting that in a homogeneous system we have u(k, k') = v l / M(k)5k,k'- We use this relation- 
ship to rewrite the homogeneous u function of Eq. ( ft.71| ) in a form suitable for generalization 
to the inhomogeneous case, 

M(k, k') = (kk'5 ky + CupS^y 1 , (3.73) 

interpreting the inversion as that of a (diagonal) matrix. In the absence of cusps, we have 
seen that the homogeneous Jastrow factor may be obtained from the inhomogeneous one by 
replacing 
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M k)k , = / e^'y^l(r)d 3 r (3.74) 

by co>p5k,k'- As we now wish to extrapolate from a homogeneous Jastrow factor to an inho- 
mogeneous one, we do the opposite and replace tu p 5k,k' = (^p^k.k') 1 ^ 2 by the matrix square 

1 /2 J-Tf-J 

root M k ' k r. Eq. ( |3.73| ) then becomes 



U (k, k') = (kk'Sw + cm^v 1 



4ttC ( r C M kk 



1/2 \ -1 

{kk>y \ 6 ^' + J • (3 - 75) 

The matrix to be inverted is no longer diagonal, but remains Hermitian and positive definite. 

If we make the reasonable assumption that the Fourier series for cUp(r) converges fairly 
rapidly, the elements of the matrix M are constant along the diagonal and fall off as we 
move away from the diagonal. For large k and k', this guarantees that tt(k, k') is dominated 
by the l/{kk') 2 ~ l//c 4 prefactor, generating a cusp. For small k and k' we have 



which is the RPA result. 

The u function of Eq. ( |3.75 ) therefore interpolates smoothly between the anisotropic 



long-range correlation term derived from the inhomogeneous RPA and the cusp at short 
range. In fc-space we now have a continuous crossover from collective behavior to two- 
particle scattering instead of a sudden and rather unphysical cutoff at k c . 

G. The one-body term 

The introduction of the cusp modifies the k < k c Fourier components of the RPA u 
function and introduces nonzero Fourier components with k > k c . In addition, it makes the 
u function spin dependent, suggesting that we need a spin-dependent one-body term. We 



therefore generalize our expression for x(k) from Eq. (|3.65|) by extending the wave vector 



sum to include components with k > k c and introducing a sum over spin indices, 
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XT (k) = -5> n (k,k')n T (k') +« n (k,k')n i (k')] , (3.77) 

k' 

with an equivalent formula for X|(k). In a spin-unpolarized system, where n^(k.') = ^j.(k') = 
-n(k'), this reduces to 

XT (k) = - E|KT(k,k') + Un (k,k')]n(k') . (3.78) 
k' z 

In the case of a homogeneous correlation term u this further reduces to 

XT (k) = ~W K T (k) + u n (k)]n(k) , (3.79) 



as first proposed by Malatesta et a/.0 

IV. RESULTS 

This section assesses the effectiveness and accuracy of QMC trial wave functions con- 
taining RPA Jastrow factors. Section |1V A] describes the systems studied and explains how 



the results are presented; Sec. |1V B| considers homogeneous systems; and Sec. |1V (J| looks at 
inhomogeneous systems. 

All the results were obtained using trial wave functions of the standard Slater- Jastrow 
form, where the spin-up and spin-down Slater determinants were constructed using accurate 
LDA orbitals. The Jastrow factor contained two- and one-body terms, u(ri,rj) and x( r i)> 
of various different types. Note that from now on we drop the w(rj, r^) self-interaction terms 
from the Jastrow factor. This is equivalent to altering x{ r i) by a negligible amount (ca. 5%). 

The energy averages discussed in the rest of this section were all accumulated using 
samples of 10000 statistically independent configurations of all the electrons; the quoted 
errors are standard deviations of the mean of the local energy. 



A. System geometry 



In one-electron bandstructure calculations the energy eigenf unctions ^k(r) may be ob- 
tained by solving the one-electron Schrodinger equation within a single unit cell subject 
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to Bloch boundary conditions. Properties of the infinite crystal may then be obtained by 
evaluating Brillouin zone integrals (i.e. averaging over boundary conditions). In QMC cal- 
culations, however, we have to take account of all the electrons simultaneously, and it is no 
longer possible to reduce the problem to a single unit cell. Instead, we consider a model solid 
consisting of a finite simulation subject to periodic (not Bloch) boundary conditions. The 
simulation cell is made as large as possible to minimize the finite-size errors, and normally 
consists of several primitive unit cells containing a few hundred electrons in total. 

The simulation cell we choose is the Wigner-Seitz cell of a face-centered cubic (FCC) 
lattice. Since we use periodic boundary conditions, any electron that moves out through 
one face of this cell immediately re-enters through the opposite face. The simulation-cell 
Hamiltonian obeys the same periodic boundary conditions, and hence a periodic model 
potential energy is required; we use the potential energy per cell of an infinite lattice of 
identical copies of the simulation cell. Since the Wigner-Seitz cell of an FCC lattice is close 
to spherical, the interactions between electrons in neighboring copies of the simulation cell 
are smaller than for most other geometries, which helps to reduce the finite-size errors. It 
is important not to confuse the simulation-cell lattice with the actual lattice structure of 
the solid; the simulation cell may contain many primitive unit cells, and these need not be 
face-centred cubic. The lattice vectors of the FCC simulation-cell lattice will be denoted 
by A 1; A 2 , and A 3 , and the corresponding body-centered cubic (BCC) reciprocal lattice 
vectors by B 1; B 2 , and B 3 . 

For reasons of computational efficiency, we have chosen to study electron gas systems 
subject to external potentials that vary along the B 3 direction only, 



where V = 1 in atomic units. Since B 3 is a reciprocal lattice vector, this choice ensures 
that the potential has the same periodicity as the simulation cell. The electron density and 
X functions, which also vary only in the B 3 direction, share this periodicity. 

The one- and two-body terms must also respect the periodic boundary conditions applied 




(4.1) 



34 



to the simulation cell. This implies that analytic Jastrow factors based on the u function 
of Eq. ( |3.69| ) must be made periodic by including contributions from all the electrons in 
a periodic lattice of identical copies of the simulation cell. Since the analytic u function 
decays like 1/r at large r, the sum of contributions is evaluated using Ewald summation 
techniques. Numerical Jastrow factors calculated from the inhomogeneous RPA are periodic 
by construction. 

The next two subsections contain a number of figures showing electron densities, x 
functions, and u functions. Charge densities and x functions are plotted along the B3 
direction from one side of the simulation cell to the other. The inhomogeneous two-body 
term u(ri,r 2 ) is more difficult to represent. We have chosen to fix the position ri of the 
first electron, while sweeping r 2 along the B 3 direction on a line passing through ri. Figure 
|l] shows the three positions of the first electron considered. Note that because A 3 and B 3 
are not the same, moving along a line parallel to B3 does not bring you back to a point 
equivalent to (i.e. differing by a lattice vector from) the starting point until you have passed 
through three layers of simulation cells. This means that although u(ri,r 2 ) always has the 
full periodicity of the simulation cell, this is not always apparent from the plots. Note also 
that all Jastrow factors are symmetric on interchange of ri and r 2 . 



B. Homogeneous systems 

The FCC simulation cell considered in this section held a uniform electron gas of 61 
up-spin electrons and 61 down-spin electrons, giving 122 electrons in total. The density 
parameter r s was equal to 2, corresponding to a Fermi wave vector /c^=0.96aQ 1 . Two 
different Jastrow factors were considered: 

a. Homogeneous RPA without cusps. The homogeneous RPA theory suggests using a 
correlation term of the form, 

^r,) = ^E^ ik(ri - rj) - ( 4 - 2 ) 

v Up k<k c K 
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We saw in Sec. |III A| that for typical metallic densities the cutoff k c is set to a value compa- 



rable to the Fermi wave vector kp. 

b. Homogeneous RPA with cusps. In Sec. [11 B| we saw how a cusp may be introduced 



by adding an exponential factor to the k c —KX> limit of the homogeneous RPA u function, 

— (l-e-^/^) , (4.3) 



u aiOj yij, 



UlpTij 



where F aiaj is chosen appropriately. 

Table | shows the local energy averages and standard deviations obtained in VQMC 
simulations using these two Jastrow factors. For comparison, we also show results obtained 
using a "Hartree-Fock" trial function including up- and down-spin Slater determinants of 
LDA orbitals (in this case plane waves) but no Jastrow factor. The introduction of an RPA 
Jastrow factor without a cusp lowers the calculated energy considerably but has little effect 
on the standard deviation. The introduction of the cusp lowers the energy greatly and also 
reduces the standard deviation. It is clear that the presence of the cusp is vital if accurate 
total energies are to be obtained. 

The cuspless RPA results shown in Table | were obtained using a value of k c equal to the 
Fermi wave vector hp, but we also investigated the limit as k c tends to infinity, in which case 
u(rij)—*—l/(u p rij). Since we know that the description of screening in terms of collective 
plasmon modes is invalid at short distances, it was no surprise that this limiting form gave 
very poor results. The residual two-electron interactions in the short-range Hamiltonian lead 
to a cusp-like behavior in the wave function at close distances, not a —l/(u p rij) divergence. 

In Sec. |riF]we saw how the u function with a cusp from Eq. ( |4.3|) may be represented 



as a Fourier series, 

" (k) = WF(FT^) • (4 ' 4) 

where C — 1 for antiparallel spins and C = 1/2 for parallel spins. We can investigate the 
usefulness of this representation by cutting off the series at a wave vector k n and varying k n 
to see how fast the calculated VQMC energy approaches the >oo limit. The hope is that 
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with a reasonably small /c-space cutoff we will be able to produce a Jastrow factor that gives 
essentially the >oo energy. Figure |2] shows the convergence of the energy graphically. It 
can be seen that a cutoff of fc n =3.95ag 1 produces a wave function with the same energy (to 
within statistical uncertainties) as an infinite cutoff. The 3.95a,7 1 cutoff is small enough to 
be computationally feasible, and so there is no difficulty in representing the cusp in fc-space. 
Note also how the standard deviation tends to zero as the energy improves. 

C. Inhomogeneous systems 

As mentioned above, the inhomogeneous systems we consider have a background po- 
tential that varies in one dimension only. The LDA electron density of the unpolarized 64 
electron simulation cell considered in this subsection is shown in Fig. [| This system is 
strongly inhomogeneous (for comparison, a typical interatomic distance in a solid is ~6a ; 
this is roughly the distance between the trough and the peak of the charge density). The 
average electron density is the same as that of a uniform system with r s =2 and Fermi wave 
vector kp=0.96aQ 1 . 

In addition to investigating the influence of the cusp, as in the homogeneous case, we 
must also now investigate the effects of the one-body x function and compare the accuracies 
of homogeneous and inhomogeneous u functions. We therefore split this section into three 
subsections: 

1. First we look at the pure (i.e. cuspless) homogeneous RPA Jastrow factor (which has no 
X function) and compare it with the pure inhomogeneous RPA Jastrow factor (which 
does have a x function). 

2. Second we investigate the effects of adding cusps to these two correlation factors. 

3. Finally we add an ad-hoc one-body x term to the homogeneous RPA Jastrow factor. 
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1. Inhomogeneous RPA without cusps 



The inhomogeneous RPA Jastrow factor considered here is the one derived in Sec. [Tl l~D| 
(Eqs. ( |3.61| ), ( |3.62| ), and (|3.63|) ), which includes both u and x functions. As always in this 
work, the matrix M is constructed using the LDA density, and the Slater determinants con- 
tain LDA orbitals. The homogeneous Jastrow factor includes the u function from Eq. ( |4.2| ) 
but no x- 111 both cases the cutoff k c is set equal to the Fermi wave vector, kp=0.96aQ \ of a 
homogeneous system with the same average electron density as the inhomogeneous system. 

Table [TJ compares the VQMC energies and standard deviations calculated using the ho- 
mogeneous and inhomogeneous Jastrow factors. It is clear that the inhomogeneous Jastrow 
factor is much the better of the two. The reason is apparent from Fig. which demonstrates 
that the inhomogeneous Jastrow factor, which has a built-in one-body term, produces a near 
optimal density. Figure |^ also shows that the homogeneous RPA Jastrow factor (which has 
no one-body term) gives a very poor electron density. This explains why the corresponding 
VQMC energy is so poor 

In Fig. [| we plot the inhomogeneous RPA two-body term. Both inhomogeneity and 
anisotropy can be seen. To aid understanding, Fig. [5] shows the Jastrow factors of three 
different homogeneous systems, the constant densities of which correspond to the local den- 
sities at the central positions of plots A, B, and C, respectively. It is clear that the three 
inhomogeneous u functions shown in Fig. ^ are much more similar than the three homo- 
geneous u functions shown in Fig. ||]. This shows that the inhomogeneous RPA u function 
is not well approximated by a local-density-like approximation based on the homogeneous 
RPA.EI In the low density region (position C), in particular, we find that the inhomogeneous 
Jastrow factor in Fig. |4] is suppressed relative to the "local density" version of Fig. [|. 

At point B, the charge density around the electron is asymmetric, and this is reflected in 
the anisotropy of the u function. The anisotropy is such that the u function is stronger on 
the side where the density is lower; this is consistent with the idea that the RPA screening 
is more effective where the electron density is high. 
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In summary, we find that the Jastrow factor has a larger range in the low density regions 
than in the high density regions, consistent with the "local density" picture of Fig. ^ and 
with the physical expectation that screening should be more effective at high densities. This 
interpretation also explains the sign of the anisotropy: the u function is weaker on the high 
density side where the screening is more effective. We find, however, that the range of 
variation of the inhomogeneous Jastrow factor is much smaller than predicted by the "local 
density" picture. This suggests that a homogeneous Jastrow factor defined by the averaged 
density of the underlying inhomogeneous system may not be such a bad approximation. 

2. Inhomogeneous RPA with cusps 

A cusp may be added to the inhomogeneous RPA u function using the Fourier-space 
method explained in Sec. [Ill F| . The results below are obtained with a Fourier cutoff k n of 
4.95<2o ; Fig. ^ suggests that this is large enough to represent the cusp accurately. Since 
we choose not to change the relationship between the u and \ functions, Eq. ( |2.15| ), the 
introduction of the cusp also modifies the one-body % function. 

The addition of the cusp to the inhomgeneous RPA Jastrow factor reduces the calculated 
VQMC energy from —13.32(4) x 1CT 2 eV per electron to —15.81(1) x 1CT 2 eV per electron. 
This is the best variational estimate of the energy we were able to obtain using any of the 
Jastrow factors considered in this paper. The addition of cusps to the homogeneous u func- 
tion does not introduce a one-body term and so the density obtained using the homogeneous 
RPA Jastrow factor is still poor. Energies calculated using the homogeneous RPA Jastrow 
factor therefore remain much worse than energies calculated using the inhomogeneous RPA 
Jastrow factor. 

Figure |6|, which is analogous to Fig. f|, shows the inhomogeneous RPA u function after the 
cusp has been added. It is clear that the addition of the cusp greatly reduces the amount of 
inhomogeneity and anisotropy. Despite the fact that the system is strongly inhomogeneous, 
the cusp acts as such a stringent constraint that «(rj, i\,-) is close to homogeneous. Although 
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the inhomogeneity derived from the RPA must persist when |r$ — r j | is large enough, the 
crossover length, 2ir/k x , corresponding to the average density r s =2, is comparable to our 
system size. This implies that the form of the u function is largely determined by the cusp 
throughout our system. If we had studied larger systems we would have seen the RPA 
reassert itself at large |rj — r^- 1 , but previous work on numerical trial function optimizationil 
and finite-size errorJiH has shown that the behavior of the u function at such large values of 
— Tj\ has very little effect on the total energy. The fact that the inhomogeneous u function 
becomes so homogeneous once the cusp has been added may explain the surprisingly good 
performance of the homogeneous u functions used in most QMC simulations of solids. 

3. Other one-body terms 



In this subsection we compare the quality of analytic one-body terms based on Eq. ( |2.15| ) 
and numerical one-body terms obtained using variance optimization. 

We have already explained that we always construct the \ function appearing in the 
inhomogeneous RPA Jastrow factor from the u function and density according to Eq. ( |2.15p . 



Although Eq. ( ^JEj) was derived within the RPA, we assume that it holds unaltered even 
after the spin-dependent cusps have been added to u. This assumption proves very successful 
in practice, yielding \ functions that are not significantly worse than those computed (at 
much greater cost) using numerical variance optimization. 

When adding a x function to the homogeneous u function of Eq. ( }4.3|) we have two 
options: we could try using Eq. ( [2.15| ) again, or we could use variance minimization. Since 
u is now a function of | r j — 1 only, Eq. ( |2.15| ) reduces to Eq. (|2.16|) , which was first derived 



by Malatesta et al. EH The impressive accuracy of Eq. ( [2.16|) is shown in Fig. [7|, where 
we compare the analytic x function with one obtained using additional numerical variance 
optimization. We see that the two x functions (and hence the two electron densities) are 
very similar. Table JTT| shows that the difference between the two VQMC energies is smaller 
than the statistical noise in our simulations. 
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Table [IV| shows that the inhomogeneous RPA (including the x function from Eq. (|2.15|) 
and an approximate Fourier representation of the cusp) yields marginally better results 
than the homogeneous RPA (including the \ function from Eq. ( |2.16| ) and an approximate 
Fourier representation of the cusp). Comparing these results to the Ewald (i.e. k n —>-oo) 
results from Table pTT| , we see that including the exact cusp reduces the variance but has 
only a negligible effect on the energy. The inhomogeneous numerical Jastrow factor yields a 
slightly better energy than the Ewald Jastrow factor. The improvement is barely statistically 
significant, but suggests that the inhomogeneity of the RPA correlations does produce a 
slight improvement. 

V. CONCLUSIONS 

Our aim was to better understand the physics underlying the Jastrow factors used in 
QMC simulations of solids, and to derive improved Jastrow factors for strongly inhomoge- 
neous systems. We began by reviewing Bohm and Pines' RPA treatment of the homogeneous 
electron gasS and generalizing it to the inhomogeneous case. The result of this analysis was 
a Slater- Jastrow trial wave function containing an anisotropic inhomogeneous Jastrow factor 
expressed as a double Fourier sum. The optimal orbitals appearing in the Slater determi- 
nants were shown to be close to Hartree-Fock or LDA orbitals, even though these theories 
do not include Jastrow factors. 

The RPA describes the long-range electronic correlations accurately, but not the 
scattering-like correlations at short distances. We saw, however, that the correct short- 
range behavior determined by the cusp conditions may easily be imposed on any Jastrow 
factor represented in fc-space. When the inhomogeneous RPA result is modified in this way, 
the result is a parameter-free Jastrow factor with the correct short and long-range behavior. 

For systems of a few hundred electrons and an external potential varying in one dimension 
only, we showed that trial functions incorporating modified RPA Jastrow factors are both 
accurate and computationally tractable. Since such Jastrow factors are parameter- free, 
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the time-consuming variance minimization procedure normally used to generate accurate 
X functions is not required. Surprisingly, however, the inhomogeneity of the two-body u 
function yields little benefit in the system we studied, producing an energy only two standard 
deviations (a barely statistically significant amount) below the best result obtained using a 
homogeneous u function. Provided we have a cusp describing the short-range interaction 
and a one-body term to mend the density, the detailed form of the u function is not very 
important. The reason is that the imposition of the cusp conditions, which fixes the gradient 
of u(ri, Yj) when |r$ — r^- 1 — > 0, washes out most of the inhomogeneity of the RPA u function 
when |rj — is small. It is the short-range correlations that have most effect on the energy, 
and so the long-range inhomogeneities that remain after the cusp conditions have been 
imposed have little effect. 

If we compare the plasmon cutoff k c (oc n 1 / 6 ) from Eq. Q3.3| ) with the wave vector 
k x (oc n 1 / 4 ) that characterizes the crossover from screening behavior to cusp-like behavior 
(see Sec. |111F| ), we see that the cusp is relatively less important in high density systems. 
This suggests that the inhomogeneity of the RPA u function may produce a more obvious 
improvement when the average electron density is both large and strongly varying. This is 
intuitively sensible, since in low density systems we expect the short-range electron-electron 
scattering described by the cusp to dominate, whereas at higher densities screening and 
collective effects should be more important. Possible candidates for high density systems 
include calculations explicitly involving the core electrons, where it is already known that the 
use of inhomogeneous u functions is advantageous.0 Other systems where one might expect 
inhomogeneities in the correlation term to become important are rare earth elements, where 
some of the valence electrons are strongly bound to the core. 

In summary, we hope that this paper has contributed a better and more general un- 
derstanding of the physics underlying the Slater- Jastrow trial functions used in most QMC 
simulations of solids. 
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APPENDIX A: PHYSICAL INTERPRETATION OF THE RANDOM-PHASE 



APPROXIMATION 

In this appendix we look at the physical interpretation of the Bohm-Pines Hamiltonian 
Hftp' after the application of the unitary transformation discussed in Sec. [Ill Q As expressed 
in Eq. ( |3.43| ) this Hamiltonian appears very complicated. However, if we define a field 

A(r) = (y)' E^ £k e ik - r (Al) 
and then calculate E(r) = — ^A(r) = -i[H£f, A(r)], we obtain: 

E(r) = -(^) 5 X>- k £ k e ikr . (A2) 
The transformed Bohm-Pines Hamiltonian may then be written in the much simpler form: 

+ ^/jE(r)]Vr + £V(r i ). (A3) 

The Hamiltonian of Eq. (|A3|) describes a set of quantum mechanical particles moving in 
a reduced external potential potential V(r) and interacting with a quantum mechanical 
longitudinal electromagnetic field A(r). The kinetic energy of the field is associated with the 
energy density of the electric field in the usual way. The interactions between the particles 
and the field are also described via the standard coupling to the momentum operators. 

The physical interpretation of the transformed subsidiary condition, f2£ ew |$ ncw ) = 0, 
where fi£ ew is as given in Eq. fl3.45| ), also becomes much clearer when re-expressed in terms 



of the new fields; it simply sets the Fourier components of 

&(r) oc divE(r) - 4n [n(r) - n(r)] (A4) 

to zero. The subsidiary condition thus ensures that the electric field is related to the density 
of the particles via Gauss's law. Crucially, we found that this constraint is automatically 
satisfied in the ground state, provided the tt° have been chosen correctly. 
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The RPA decouples the electronic and field variables by neglecting the p • A terms and 
treating the quadratic A 2 terms only approximately by averaging the electron positions. 
The field part of the decoupled Hamiltonian is then harmonic: 





A( r y 







n(r)d 3 r H 



I 


E( r y 







d 3 r . 



(A5) 



The resulting wave function no longer obeys Gauss's law automatically, but a projection 
operator can be applied to produce a wave function that does. 

Eq. ( |A5|) looks disconcertingly simple: it appears that we have three independent har- 
monic oscillators at every point r (one each for the x, y and z components of the vector 
potential). The apparent simplicity is misleading, however, since the condition that the 
fields be longitudinal, curlA = 0, couples the oscillators at different points in space and 
reduces the number of degrees of freedom at any point from three to effectively just one. 
Dropping this condition and retaining only one scalar oscillator is equivalent to dropping 
the £k ■ £k' term in Eq. ( |3.50|) and leads to a local- density-like version of the RPA. 
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FIGURES 



Position of Electron 1 in Plot A 



in Plot C 




FIG. 1. In figures showing two-body terms, plots labeled A show u(ri,r2) as a function of T2 
for i"i fixed at the peak of the electron density; plots labeled B show u(ri,r2) for ri fixed at the 
average of the electron density; and plots labeled C show u(ri,r2) for ri fixed at the minimum of 
the electron density. In all cases, T2 is swept along the B3 direction on a line passing through n. 
The relative coordinate z measures the distance between the two electrons and thus equals zero 
when the two electrons are at the same place, irrespective of the fixed position of r\. 
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FIG. 2. The convergence of the VQMC energy as a function of the cutoff k n used in the 
truncated Fourier series representation of the u function from Eq. ( |4.3| ). The results are for the 
uniform system considered in Sec. IV B, for which fop = 0.96CIQ 1 . The dotted line shows the 
calculated value of the energy when k n =oo (the standard deviation of this result is too small to 
show here). 
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12.0 



FIG. 3. The electron density of the strongly inhomogeneous 64 electron system considered in 
Sec. [V C[ The LDA density (solid line) is compared to the densities obtained using the homoge- 
neous RPA (dotted line) and the inhomogeneous RPA (dashed line); the z-axis lies along the B3 
direction. 
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FIG. 4. The inhomogeneous RPA u function with no cusp for three different positions of 



the fixed electron. The results are for the inhomogeneous system considered in Sec. IV C. The 
definition of z and the positions of A, B, and C are explained in Fig. |l]. The Jastrow factor is 
stronger in the low density region. 
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FIG. 5. The RPA u functions for three different uniform electron gases, the densities of which 
are equal to the densities at points A, B, and C of the strongly inhomogeneous 64 electron system 



considered in Sec. IV C. The definition of z and the positions of A, B, and C are explained in 
Fig. |l[ The homogeneous u functions are of course isotropic. 
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FIG. 6. The inhomogeneous RPA u function with a cusp for three different positions of the 
fixed electron. The results are for the strongly inhomogeneous 64 electron system considered in 
Sec. [TV Q . The definition of z and the positions of A, B, and C are explained in Fig. [I]. The 
addition of the cusp has much reduced the inhomogeneity and anisotropy observed in Fig. |4|. 
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FIG. 7. Comparison of the x function (solid line) obtained from Eq. ( |2.16| ) with one (dashed 
line) obtained using an additional variance minimization. The results are for the strongly inho- 
mogeneous 64 electron system considered in Sec. IV Q , using the homogeneous Ewald-summed 



Jastrow factor with cusp. The corresponding energies, shown in Table |ni| are equal to within the 
statistical error. 
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TABLES 

TABLE I. VQMC local energy averages and standard deviations of the uniform system consid- 



ered in Sec. IVB , Results for three different trial wave functions are shown. The HF trial function 
has no Jastrow factor. The RPA results use the "pure" RPA u function from Eq. ([4.2|). The best 
energies are obtained using the RPA u function with a cusp from Eq. 





HF 


RPA 


RPA+CUSP 


Energy (10 2 eV per electron) 


2.95 


1.94 


-0.59 


Std. dev. (10 -4 eV per electron) 


2.79 


2.60 


0.50 



TABLE II. VQMC local energy averages and standard deviations of the inhomogeneous system 



considered in Sec. IV C. The HF trial function has no Jastrow factor. The RPAI trial function 
includes a Jastrow factor containing the homogeneous RPA u function from Eq. fl4.2|) but no \ 
function. The RPAII trial function uses the full inhomogeneous RPA Jastrow factor derived in 





HF 


RPAI 


RPAII 


Energy (10 2 eV per electron) 


-12.24 


-5.8 


-13.32 


Std. Dev. (10 _4 eV per electron) 


4.0 


4.4 


3.6 



TABLE III. VQMC local energy averages and standard deviations of the system considered 



in Sec. IV C. Results for three different trial wave functions are shown. The HF trial function has 
no Jastrow factor. The EW trial function includes the Ewald-summed u function from Eq. (| 
and an analytic \ function calculated using Eq. ( 2.16j ). The VM trial function uses the same 



u 



function but optimizes the Fourier components of x using variance minimization. 







HF 


EW 


VM 


Energy 


(10 2 eV per electron) 


-12.244 


-15.786 


-15.790 


Std. Dev. 


(10 _5 eV per electron) 


40.0 


8.2 


8.4 



52 



TABLE IV. VQMC local energy averages and standard deviations of the system considered 



in Sec. IV C . Results for three different trial wave functions are shown. The HF trial function 
has no Jastrow factor. The RPAI trial function includes the homogeneous RPA u function with 
Fourier components given by Eq. (4.4) plus a \ function generated using Eq. fl2l6|) . The RPAII 



trial function uses the inhomogeneous RPA Jastrow factor from Sec. HID, to which cusps have 
been added as explained in Sec. IIIF, The Fourier cutoff k n was set to 4.95oq 1 in both cases. 







HF 


RPAI 


RPAII 


Energy (10 


~ 2 eV per electron) 


-12.244 


-15.785 


-15.812 


Std. Dev. (10- 


_5 eV per electron) 


40.0 


12.2 


11.4 



53 



